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Abstract. 

The occurrence of a bounce in FRW cosmology requires modifications of general relativity. 
An example of such a modification is the recently proposed Hofava-Lifshitz theory of gravity, 
which includes a "dark radiation" term with a negative coefficient in the analog of the Fricdmann 
equation. This letter describes a phase space analysis of models of this sort with the aim of 
determining to what extent bouncing solutions can occur. It is found that they are is possible, 
but not generic in models under consideration. Apart from previously known bouncing solutions 
some new ones are also described. Other interesting solutions found include ones which describe 
a novel sort of quasi stationary, oscillating universes. 
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1. Introduction 

The standard ACDM model has solved many issues in cosmology. However, in spite of all this 
success, it also leaves a number of issues unaddressed. Perhaps the most significant is the problem 
of initial singularity, where general relativity breaks down. There have been many attempts 
to modify Einstein's theory to avoid this singularity. Some are made at classical level, some 
involve quantum effects. Examples include the ekpyrotic/cyclic model ( [H [2J, [3] ) and loop quantum 
cosmology (jll|5j[6]), which replace the Big Bang with a Big Bounce. Attempts to address these 
issues at the classical level include braneworld scenarios (0 [5]), where the universe goes from 
an era of accelerated collapse to an expanding era without any divergences or singular behavior. 
There are also higher order gravitational theories and theories with scalar fields (see for a 
review of bouncing cosmologies) . It is however fair to say that the issue of the initial singularity 
remains one of the key questions of early Universe cosmology, and the idea that it is avoided due 
to a bounce is remains an elusive (and promising) notion. As discussed below, it is clear that close 
to the singularity the Friedmann equation has to be modified for a bounce to be possible. 

In recent months much effort has been devoted to studies of a proposal for a UV complete 
theory of gravity due to Hof ava (TUJ [JTJ E3] ■ Because in the UV the theory possesses a fixed 
point with an anisotropic, Lifshitz scaling between time and space, this theory is referred to as 
Hofava-Lifshitz gravity. Soon after this theory was proposed many specific solutions of this theory 
have been found, including cosmological ones ([TJl HSl HE1 H3 [IH1 UHl 1201 EI])- It was also realized 
that the analog of the Friedmann equation in HL gravity contains a term which scales in the same 
way as dark radiation in braneworld scenarios [H] E3 US] and gives a negative contribution to 
the energy density. Thus, at least in principle it is possible to obtain non-singular cosmological 
evolution within Hof ava theory, as it was pointed out in [H] [TBI EZ] • 

Although there is presently much discussion of possible problems and instabilities of Hofava- 
Lifshitz gravity [531 [Ml 123 US] j it is still very interesting to perform a detailed investigation of 
the influence of the additional terms in the Friedmann equation in HL gravity on the existence 
and stability of a cosmological bounce. In this paper we are going to analyze how these terms 
affect the dynamics of the system using phase portrait techniques described in 32, 33 , and then 
compare the results with those valid in standard cosmology. For purpose of illustration we will 
assume that matter in the pre-bounce epoch is described by a scalar field if with a quadratic 
potential. In order to concentrate on modifications created by the "dark radiation" terms, we will 
set cosmological constant A = 0. Such scenario may be also considered as an approximation to 
a general case with non-vanishing A, valid in the regime of small scale factor a, when standard 
curvature and cosmological constant terms (vanishing for A = 0) become negligible. Thus, the 
present analysis can be regarded as an exploration of the cosmologies with modified equations 
of motion, where the particular modifications considered are inspired by Hofava cosmology. In 
particular, the main question addressed is the impact of these modifications on the existence of a 
bounce, which is otherwise not possible. 

Related analysis of Hofava-Lifshitz cosmology have recently appeared in [27] and [28] , which 
we become aware of while this work was being typed. Those papers address a somewhat different 
set of issues from what we have pursued. The analysis presented in [27) and [28] consider the full 
4-dimensional phase space of HL cosmology. The results presented here focus on the region close 
to where the scale factor vanishes, which admits a critical simplification: the number of dynamical 
equations under study can be reduced from 4 to 3 (as discussed in more detail in Section 3). This 
makes it possible to visualize the possible phase space trajectories in a 3-dimensional space. 

The structure of this note is following: in Section 2. we briefly sketch Hofava-Lifshitz gravity 
and cosmology. In Section 3. the possibility of bounce is discussed. In Section 4 we discuss phase 
portraits of the discussed system of equations and describe different families of phase trajectories. 
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2. Hofava-Lifshitz cosmology 

The metric of Hofava-Lifshitz theory, due to anisotropy in UV, is written in the (3 + l)-dimensional 
ADM formalism: 

ds 2 = -N 2 dt 2 + g tj {dx l - N i dt)(dx j - N 3 dt), (1) 
where N, iV$ and g^ are dynamical variables. The action of Hofava-Lifshitz theory is Q2] 

/ dtd 3 x{C + &), (2) 



I 



r 17,n\ 2 (k i «V(Ai?-3A 2 ) 

^ = JgN^iKvK'- XK ) + g(j - ^ 



A-v^j 32(1 -3A) * "2^ ZljZ 

where ify = ^ [jSij — ^(i-^Yj)] ^ s extrinsic curvature of a spacelike hypersurface with a fixed time, 
a dot denotes a derivative with respect to t and covariant derivatives defined with respect to the 
spatial metric gij. Moreover 

flUJ 2 

Zij = Cij ^-Rij- (3) 
k 2 , A, /i, oj and A are constant parameters and the Cotton tensor, C^ , is defined by 

C ij = £j fei Vfc ^ _ l^i^J = c <*iv fc i^, - ^e lfcj a fc i?. (4) 

In ([2]), £ is the kinetic part of the action, while C\ gives the potential of the theory in the so-called 
"detailed-balance" form. 

Matter may be added by introducing a scalar field (p ([21 US]) with energy density p and 
pressure p. The action for matter is 

I m = J dtd 3 x^gN£ m . (5) 

The matter Lagrangian C m depends on the scalar matter field ip and the 4-dimensional metric. 
r ■ 

C m = (^(^ - Nid w) - V (V)) (6) 

This allows to define the energy density and pressure of the scalar field in the following way: 

p = 3 -^lp 2 + V(p), (7) 

p = - V(p) (8) 

In numerical calculations presented further on a specific form of the scalar potential will be assumed 
(see eqn. (14l). 

Comparing the action of Hofava-Lifshitz theory to the Einstein-Hilbert action of general 
relativity, one can see that the speed of light, Newton's constant and the cosmological constant 
are 



k 2 fi r^r~ „ k 2 c 3«v 2 a 2 



4 V 1-3A' ° 32^' Ae 3A-1 32' (9) 
respectively. Setting dynamical constant A = 1, reduces the first three terms in ^ to the usual 
ones of Einstein's relativity and matter Langrangian in ([5| to the usual scalar field action in curved 
space-time. 

The equations for Hofava-Lifshitz cosmology are obtained by imposing condition of 
homogeneity and isotropy of the metric. Precisely, the equations of motion are obtained by varying 
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the action ^ with respect to N, a, and ip, and setting N = 1 at the end of the calculation, leading 
to 

2 K 2 p ft 4 // 2 A k K 4 fi 2 / 2 k 2 \ 

= 6(3A - 1) + 8(3A - l) 2 a 2 _ i6( 3A - 1)2 ^ + ) ' ( 10 ) 

4(3A-1) 8(3A-l) 2 a 2 32(3A-l) 2 a 4 ' 1 J 

and also equation of motion for the scalar field: 

+ 3^ + ^-^ = 0, (12) 

where H = a/a, a prime denotes the derivative with respect to scalar field if. The significant 



new terms in the above equations of motion are the (1/a )-terms on the right-hand sides of (10 1 



and (11). They are reminiscent of the dark radiation term in braneworld cosmology [29] and are 
present only if the spatial curvature of the metric is non- vanishing. 

Equations ( |To| ) and ( |TT] ) show different behavior for different ranges of the A-parameter: 
A>l/3,A = l/3 and A < 1/3. It was shown in [30 that 1 > A > 1/3 leads to ghost instabilities 
in the IR limit of the theory. Solution to this problem proposed in [31 results in instabilities 
re-emerging at UV region. However, this range of A is exactly the flow-interval between the UV 
and IR regimes. Thus the only physically interesting case that remains, allowing for a possible 
flow towards GR - at A = 1 - is the regime A > 1. Region A < 1/3 is disconnected from A = 1 
and cannot be included in realistic considerations. 

Thus we will remain in the phenomenologically range relevant 00 > A > 1. In this case the 
value of A - a dimensionless coupling of the theory - may be included in rescalling of the parameter 
k. In general, A runs - logarithmically in the UV - and may eventually reach one the three IR fixed 
points (pj]): A = 1/3, A = 1 or A = 00, the first one excluded by the existence of instabilities. 

3. Existence of bounce 



New terms in the cosmological equations introduce the possibility of a bounce. The form of (10), 
with k = ±1 implies that it is possible that H = at some moment of time. This is a necessary 
condition for the realization of the bounce. It was pointed out in [14] . that it may happen in the 
presence of matter, at the critical time £*, a = a*, when the critical energy density is equal to 

3k 2 u 2 / A k A 2 lk 2 \ 
P = P* = wv^7 + ■ (13) 



3A - 1 V 4 a 2 8 8 a 4 , 
which is determined by the couplings of the theory. 

From the continuity equation it follows that at the bounce point H > 0. Therefore a transition 
from a contracting to an expanding phase may be possible. It was shown in 22 that the necessary 
condition for a cosmological bounce is that the energy density of regular matter increases less fast 
than a~ 4 as the scale factor decreases and — p) > . 

We begin our considerations during a contracting phase. At the beginning the scale factor is 
quite large and the contribution of dark radiation to the total energy density is quite small. As 
the universe contracts, the energy density increases and the scale factor decreases rapidly. When 
a critical density is achieved, a big bounce is about to take place. 



One would expect that near the bounce the leading term in ( 10 ) and ( 11 ) would be the dark 
radiation one, with curvature and cosmological constant terms neglectable small. Actually the 
latter terms generally vanish when HL cosmological constant A = 0. In case of non-vanishing 
A these two terms may be neglected when the scale factor a is sufficiently small. Specifically, 
assuming for a moment an equation of state of the form p — wp with constant w, it is well known 
that H 2 , H and p scale as a~^ 1+w \ Therefore we may keep the density term and omit the 
curvature term ~ 1/a 2 if w > — |. In the case of a quadratic potential considered below (for 
which w ^ const, so the above argument does not directly apply) we have checked numerically 
that in all bounce scenarios discussed in this paper, this approximation is valid near bounce point 
(up to lO" 7 ). 
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We will model the matter sector in this pre-bounce epoch by assuming it is described by a 
scalar field (p with a potential 

V(<p) = \m\\ (14) 

For calculational simplicity we put m = 1. 

This way and inserting a — 2/k 2 , we have the following equations modeling bounce in the 
Hofava-Lifshitz cosmology. 

The value of k 2 may be expressed in terms of cosmological constants |9]), we will work in units 
such that 8ttG — 1 and c = 1. Then 

k 2 = 32ttGc, (17) 

and the values of ^ are left arbitrary. Therefore the Friedmann equations take the following form 
near the bounce: 

1 , 2u 2 k 2 , , 

" = - S^T* + (3A^W < 18) 

Additionally, completing dynamics of the system, there is the equation of motion for the scalar 
field and the definition of the Hubble parameter: 

(p = - - tyH, (20) 

a = aH. (21) 

The value of the parameter /i may be kept arbitrary. This parameter does not alter solutions of 
the system ( 18][21 ), but it specifies values of a on an obtained trajectory. 



4. Phase portrait 

Phase space 

The local geometry of the phase portrait is characterized by the nature and position of its critical 
points. These points are locations where the derivatives of all the dynamic variables, i.e. the 
r.h.s. of ( 33p5 ), vanish. Moreover, they are the only points where phase trajectories may start, 



end, or intersect. They can also begin or end in infinity, and then - after a suitable coordinate 
transformation projecting the complete phase space onto a compact region (so called Poincare 
projection) - there may be well defined infinite critical points. The set of finite and infinite critical 
points and their characteristic, given by the properties of the Jacobian matrix of the linearized 
equations at those points, provides a qualitative description of the given dynamical system. 
Dynamics of our system is described by the following set of first order ODE's: 

u = (p (22) 
u = - — ^ - SuH, (23) 
a =aH (24) 



plus the constraint equation 

1 

3(3T-1) V ~ ' ^ ; (3A-l) 2 a 4 



H 2 = TTToT" tt( u ' 2 + < /' 2 ) - (26) 
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If spatial curvature k — one may consider a 2-dimensional subsystem: 

u = ip 

u = — — -<p — 3uH 7 

2 



H 



3A - 1 
1 



3A - 1 

with a constraint equation 



- u 



(27) 
(28) 

(29) 
(30) 

3( 3A-ir - (31) 

If k 7^ one may also consider a subsystem on variables (ip,u,H), obtained via reduction of the 
original system with respect to constraint ( 26 ) . Namely substituting 

M 2 fc 2 1 



H 2 



1 



(3A-l) 2 a 4 3(3A-1 



-(u 2 + ^ 2 )-h 2 



(32) 



into the equation for H and omitting equation on dynamics of a leads the following set of equations: 
u = ip (33) 

(34) 



H = 



3A - 1 
2 



ip — SuH, 



- ^) - 2H 2 



3(3A-1)^ 2' — ' (35) 
This is a reduced 3-dimensional subset of ( 22p6 ) on variables ((p,u,H). If one wants to obtain also 



dynamics of a, he needs to add to this system equation a = aH and also the constraint equation 



(26). 



In subsequent considerations we shall focus on a case k ^ (when HL corrections play 
a significant role) and a phase portrait of solution of the system ( 33p5 ) in space of (ip,u,H), 
following similar a similar procedure as that described in ( |32[ 133] ) . Reducing dimensionality of 
phase space enables 3D phase portrait visualizations. Moreover, we will discuss shortly also a case 
k = 0, which play a role of a limiting case of k ^ dynamics. 

We start by rewriting equations ( 33|[35 1 in terms of the variables 



x = ip; y = (p; z = -, 
a 

which gives three "evolution" equations 

x = y, 

9 

3yz, 



3A- 1 
2 



-M 2 



V 



)-2z 2 



(36) 

(37) 
(38) 

(39) 



3(3A-1) V ~ 2 

The space of solution of the above dynamical system is a 3D region of the phase space 
(x,y,z). This region is bounded by a 2D surface defined by a constraint equation (31) - space 
of trajectories of a flat universe (k = 0). This limiting surface is a double cone, with the upper 
branch corresponding to expansion and lower one to contraction. Those two branches connect at 
a point: (0,0,0), which is a critical point (see below). Hence there are no trajectories passing from 
one branch of the cone to the other. For k — ±1 all trajectories lie between the branches of this 
cone. This cone is also a limiting surface for trajectories with large a. The further a trajectory 
lies from this cone, the smaller are the values of a along it. 

Moreover, with the value of A varying between 1 and oo, the double cone reduces to the 
surface z = when A — oo. Thus for finite values of A the whole dynamics of the system is 



contained within the double cone z 2 — ^ {x 2 + y 2 ), and for A = oo all phase point lay on 
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the surface z = 0. In the latter case the phase dynamics is flat, the Universe is static, and there 
are no bounce points. Such a situation is not interesting when searching for a bounce. For finite 
A the qualitative description of the system does not depend on the specific value of this coupling 
constant, as it corresponds to the angle of the limiting cone and "shrinking" trajectories between 
double cone's branches. Therefore in order to simplify further calculations we will set the value 
A = 1. 

The bounce happens when a phase trajectory passes from the region z < to region z > 0, 
intersecting the plane z — 0. At the crossing point z > must hold. Equation [39] implies that 
this happens if the crossing point is contained between lines y — y/2x and y = — \/2x laying on 
the plane 2 = 0. Those lines are the z — section of an elliptic cone g (x 2 — — z 2 = 0, whose 



interior consists of trajectories with i > (eq. (39)). The area outside this cone is filled with 
trajectories along which z < 0. 

To find the finite critical points we set all right-hand-sides of equations ( 37p9 ) to zero. This 
gives rise to the conditions 

x = y = z = 0. (40) 

Stability properties of this point are determined by the eigenvalues of the Jacobian of the 
system ( 37p9 ). More precisely, one has to linearize transformed equations ( 37p9 1 at each point. 



Inserting x — xq + 5x, where x = (x,y,z), and keeping terms up to 1st order in Sx leads to an 
evolution equation of the form 8x — ASx. Eigenvalues of A describe stability properties at the 
given point. 

At the finite critical point O = (0,0,0), the matrix A has 2 purely imaginary eigenvalues, 
which implies there are closed orbits in the xy-plane encircling the z-axis, i.e. point O lays on 
a center line surrounded by closed orbits. It is interesting to note, that in general case A ^ 1 



eigenvalues of A are following: (y 1 _ 3A , — y 1 _ 3A , 0). Hence for A > 1/3 point O lays on a center 

line, for A < 1/3 it is a saddle. Nonetheless, the latter case is excluded form our considerations, 
as we explained in the introduction. 

To find critical points that occur at infinite values of the parameters we rescale the infinite 
space (x, y, z) into a finite Poincare sphere by means of the coordinate change: 

X (41) 

(42) 

(43) 



1-r' 
Y 

z 
T^7' 



where 



X = r sin cos ^, (44) 

Y = r sin 9 sin ip , (45) 

Z = r cos 6, (46) 

r 2 =X 2 + Y 2 + Z 2 . (47) 

We shall use both Cartesian coordinates (X, Y, Z) and spherical ones: (r,6,ip). We also rescale 
the time parameter t by defining new time parameter T such that: dT = dt/(l — r). In these 
coordinates our phase space is contained within a sphere of radius one - infinity corresponds to 
r = 1. 

This is a conformal transformation, hence the limiting cone for phase trajectories is Z 2 = 
\{X 2 + Y 2 )-. all physical trajectories are contained within this cone. Bounce points are located 
on the plane Z = within the region bounded by lines Y = V%X and Y = -V2X. The 
region containing trajectories with Z > (i.e. with H increasing) is bounded by an elliptic cone 

Ux 2 -^)-z 2 = o. 
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Table 1. The properties of the infinite critical points. 



After Poincare transformation, equations ( 37][39 ) take the following form, written in the 
spherical coordinates (r, 0, ip): 

r = ^ ~ 1 ^ cos 9 [82 + 14 cos 29 - 42 cos 26 
48 

+21 cos 2(0-0) +21 cos 2(0 + 0)], (48) 
6' = ^rsin6»(5 + 7cos26»)(l + 3cos20), (49) 
6' = r — 1 — 3r cos 8 cos 6 sin 6. (50) 

The form of the above equations is similar to the ones obtained in [32l [33] . Taking limit r = 1 
and putting r.h.s. of equations for 9' and 0' to zero, we find 12 solutions for 8, 6 at the Poincare 
sphere, shown in the Table 1. As we can see, there are 4 saddle points (more precisely saddle lines 
with end points at Si, S%, S3, S4). In the contracting part of the phase portrait (z < 0) there are 
two attracting nodes A 2 and A4 and two repulsing lines starting at R\ and R3. Hyperbolic areas 
near the nodes are bounded by repulsing lines which play role of separatrices. The expanding part 
is a mirror ("reversed in time") of the contracting one. Stability properties of those point does 
not depend on the value of parameter A, unless it is in the range (1/3, 00). 

Stability properties of infinite critical points are described in the Table 1, their position in 3D 
phase space, on a Poincare sphere, is shown in Fig. 1. 



4-2. Trajectories 

When spatial curvature k — then phase trajectories lay on the limiting cone z 2 — \{X 2 + Y 2 ), 
as shown in the Figure 2. In the contracting part all trajectories start winding around z-axis, 
then some of them end at attracting node A2, some at A4. Those two families are separated by 
repelling lines with end-points at R\ or i?3, acting as separatrices. Expanding part is a mirror 
reflection with time reversed of the contracting part. 

Trajectories of non-flat universes lay inside the limiting cone of flat space. In the contracting 
part of the diagram [Z < 0), trajectories start spiraling outside from circles around the Z-axis. 
There are two families of such trajectories, separated by repelling lines ending at R\ and -R3. In 
each family, there are two possible scenarios for subsequent evolution. The first one, shown in 
Figure 3a, is to end at an attractor node {A2 or A4), which also lays in the contracting part of 
the phase diagram. On the way between O and A 2 or A^, a trajectory may go up through Z = 
surface, undergoing a bounce there, and then recollapse, crossing the Z = plane again, or go 
straight to the attractor node, without bounce. In either case, the end is a Big Crunch. 
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Figure 1. Infinite critical points located on a Poincare sphere 




Figure 2. Phase trajectories for flat HL universe 



The second scenario is shown on the Figure 3c. Here, after some oscillations and H decreasing, 
the trajectories reach an attractor - a repelling line (that ends either at R\ or R3), along which they 
move until H = 0. Then they rapidly go up, crossing the Z — (i.e. H = 0) plane, undergoing 
a bounce. After that, and after a period of accelerated expansion, they reach another attractor - 
an attracting line laying in the expanding part, with endpoint at cither A± or A3. Along this line 
trajectories approach the Z-axis, winding around it. A subcase of this scenario is shown on the 
Figure 3e, where trajectories do not go through accelerated contraction and expansion, but cross 
the Z = surface during oscillations around the Z-axis. This is in fact the scenario described in 

Trajectories may also start at repelling nodes R2 or i?4 in the expanding part of the diagram. 
At those points H = 00, i.e. there is a Bing Bang there. After that and a period of extreme, 
but with decreasing rate, expansion, there are again two possible scenarios. One is shown on 
the Figure 3b, where trajectories reach an attracting line and end up winding around the Z-axis. 
Before that, some of them collapse, crossing the Z — plane, then slow down and finally stop 
contraction culminating in a bounce. Others show only slow expansion, without crossing Z = 0. 
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The last scenario is shown on the Figure 3d. Trajectories start at Big Bang points R2 or 
i?4, and after a period of slowing down expansion, reach a turning point and start accelerated 
contraction, ending at Big Crunch points A2 or A4. 

For better visualization we have gathered some described families of trajectories Figure 3f. 

Special attention has to be paid to circular motion around 2-axis. As stability properties 
of the point (0,0,0) and constraint equation (261 suggest, there may exist closed circular orbits 
laying on a Z = plane (H — 0). But they cannot. Equation (39) does not allow for this, as 

Z = is fulfilled only on a class of curves laying on X 2 — ^ 6Z 2 = 0, i.e. on a surface of 

the elliptic cone mentioned before. Yet numerical simulations exhibit oscillating solutions, such as 
X 2 + Y 2 = const, and Z oscillated around zero, Z > between lines Y — y/2X and Y = —y2X, 
Z < outside this region. Such a trajectory resembles deformed circle. These solutions appear 
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for sufficiently small X and Y, for larger values of X and Y numerical simulations show slow 
decreasing of the radius of this " circle" . 

Finally, note that except for the special solution shown in Figure 3e - which is the bounce 
described by Brandenberger [55], there are also other types of bounces. One, probably the most 
interesting, is shown in Figure 3c. Here a big existing universe slowly starts to contract, but later 
on the contraction becomes exponential, until a bounce is reached and an exponential expansion 
begins, which finally slows down. Another type of bounce, shown on the Figure 3a, happens again 
when a big universe slowly starts contracting, stops and goes through an expanding phase for a 
while, then recollapses and ends at Big Crunch. The last one, shown on the Figure 3b, happens 
during a transition from a Big Bang to a quasi stationary final stage (with H slowly decreasing), 
however with a bounce on the way. 

Trajectories shown on the figures discussed above are numerical solutions of the equations 
( 48][50 ). To find different bounce scenarios we investigated initial conditions: 9 — ir/2 ± 0.01, 
cj) = i|j (i= 1. . .20), each for r = j/10 (j = 1 . . . 9) and r = 0,9 + j -0.01; time in range [-20,20]. 
This procedure picked up the classes of trajectories discussed above. In general it may not be 
exhaustive in the sense that qualitatively different behavior of solutions may be possible. However 
it is sufficient for the purpose of understanding how bouncing scenarios emerge here due to the 
specific modification of general relativity which appears in Hofava's theory. 



5. Discussion and conclusions 



In this paper we have investigated the cosmological bounce in Hofava-Lifshitz gravity. Using a 
3D flow visualization technique we have found that phase portraits in the considered theory have 
a different structure than in standard cosmology. Comparing to results from the paper j33j . we 
can see that here are additional repellers (i?i and R2) in the contracting part of a phase space, 
and mirror attractors in the expanding part. Their presence allows the existence of a bounce, 
because now there are possible new families of trajectories, starting at additional repellers in 
the contracting part, and possibly ending at new attractors in expanding part, or surrounding 
the (0,0,0) point, which is now a center, compared to saddle in standard cosmology. Those are 
realizations of the bounce. One of them is the solution with oscillatory behavior described in 
( 22J); there are however additional possibilities. The most interesting one contains a period of 
rapid contraction, and - after a bounce - a period of rapid expansion, what may fit inflationary 
scenario. 

Nevertheless there are still initial conditions which lead to a Big Crunch, as shown on the 
Figures 3a and 3d, or which start at initial singularity (Fig. 3b and 3d). Hence the existence of a 
bounce is not generic for Hofava theory and depends on initial conditions. 

Another interesting class of solutions consists of quasi stationary universes. These solutions 
are described in phase space by closed orbits, winding around the critical point (0, 0, 0) - a center. 
All trajectories in the neighborhood of this point end up as closed orbits, "deformed circles". 
Equations of motion do not allow for closed orbits laying on Z = const, plane, resulting in slight 
deformation of the circular orbits. The values of H oscillate around stationary stage, for sufficiently 
small values of tp and tp. Values of scale parameter a during this evolution are much bigger than the 
regime for which our simplifications are valid. Therefore this behavior is not a feature of Hofava- 
Lifshitz theory, but of cosmologies with modified equations of motion, i.e. with the additional 
term ~ 1/a 4 in the Friedmann equations. Still, presence of this term leads to a different solution 
than induced by a negative potential as in [33] , due to different stability properties of finite critical 
points there. 

The visualizations described in this paper describe the dynamics of Hofava-Lifshitz universe 
with vanishing cosmological constant A, or HL universe with non-zero Lambda in the region of 
small scale factor a. Even in such slightly limited framework they answer the question of possible 
scenarios realizing a bounce, and whether it is generic for the theory or not. It appears not, as we 
have found solutions leading to a Big Crunch, or starting at Big Bounce, both staying within the 
regime of small a. There is also an interesting possibility of quasi stationary, oscillating universe, 
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existence of which is clearly implied by dark radiation term in Friedmann equations. 

Finally, it is worth stressing that the analysis presented here should be applicable to other 
theories which lead to modifications of the Friedmann equations. 
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